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INVERSION  OF  ISOPARAMETRIC  MAPPINGS  AND  APPLICATIONS 
OF  COMPUTER  GRAPHICS  TO  THREE  DIMENSIONAL 
FINITE  ELEMENT  ANALYSES 

ABSTRACT 

Finite  element  analyses  of  two  and  three  dimensional  structures 
often  require  the  use  of  curved  isoparametric  elements.  The  overall 
effectiveness  of  such  analyses  by  large  scale  finite  element  com- 
puter programs  is  often  handicapped  by  the  costly  and  tedious  Job 
of  verifying  input  data  and  the  lack  of  a comprehensive  graphical 
presentation  of  the  output.  This  contract  work  involved  mathematical 

research  concerning  various  approaches  to  improving  this  pre-  and 
post  processing  of  finite  element  data  utilizing  interactive  computer 
graphics.  A key  aspect  of  this  work  was  the  characterization  of 
the  existence  of  the  inverses  of  isoparametric  mappings,  as  well 
as  the  development  of  efficient  algorithms  for  the  numerical  inver- 
sion of  such  mappings. 
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INVERSION  OF  ISOPARAMETRIC  MAPPINGS  AND  APPLICATIONS 
OF  COMPUTER  GRAPHICS  TO  THREE  DIMENSIONAL 
FINITE  ELEMENT  ANALYSES 

1.  INTRODUCTION 

Finite  element  analyses  of  three  dimensional  structures  with 
curved  bounding  surfaces  most  often  require  the  use  of  curved 
"isoparametric  brick"  elements,  C63.  One  such  brick  is  illustrated 
in  Figure  1. 


r 


FIGURE  1.  The  20-node  isoparametric  brick  element 
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The  trial  functions  In  the  finite  element  (or  Ritz-Galerkin) 
method  for  such  elements  are  expressed  in  terms  of  the  curvilinear 
(or  natural)  coordinate  system  (r,s,t)  and  the  computations  of 
stiffness  matrices,  etc.  are  carried  out  in  this  same  local 
system;  not  the  global  (x,y,z)  system.  Since  the  unknowns  are, 
say,  the  displacements  at  the  nodal  points,  there  is  no  problem  in 
going  from  the  local  (r,s,t)  system  to  the  global  (x,y,z)  system. 

That  is,  built  into  the  transformation  T is  an  explicit  1-1  corres- 
pondence between  the  20  nodes  of  E and  the  20  nodes  of  the  "standard" 
brick  J : (0,1)  * (0,1)  * (0,1).  For  example,  the  geometric  des- 
cription of  £ includes  global  coordinates  (x1,y1,z1)  of  node  1 and 
? associates  node  1 with  (0,0,1)  in  the  (r,s,t)  local  system. 

However,  such  an  explicit  relation  for  points  other  than  nodes 
does  not  exist  in  general. 

This  leads  to  the  investigation  of  the  question  of  existence 
of  an  inverse  mapping  T_1:  E •*  >/  as  well  as  efficient  algorithms 
for  its  evaluation.  An  understanding  of  the  dependence  of  the 
invertibility  of  T on  the  choice  of  nodal  topology  is  imperative 
for  the  successful  construction  of  finite  element  idealization  of 
a given  structure. 

The  question  of  the  existence  of  an  inverse  itself  is  of  interest 
since  it  is  not  guaranteed  a priori.  Heuristic  techniques  for 
checking  the  existence  of  an  Inverse  which  rely  on  graphics  are 
presented  in  C 33 . 

The  overall  effectiveness  of  solving  three  dimensional 
problems  by  multipurpose  finite  element  computer  programs  is 
handicapped  by: 


-1 
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(i)  the  costly  and  tedious  Job  of  setting  up  and  verifying 
Input  data  (esp.  geometric  data),  Cl-,2,3,53, 

(11)  the  lack  of  a comprehensive  graphical  presentation  of 
the  computed  stress  distributions  for  three  dimensional  structures, 
C5,73. 

By  utilizing  newly  developed  tools  of  numerical  analysis  and 
computer  graphics,  numerical  methods  and  associated  computer  programs 
were  developed  which  can  be  used  in: 

(1)  a pre-processor  mode  to  aid  in  more  efficiently  describing 
and  verifying  the  geometric  .ir\put  and, 

(ii)  a post-processor  mode  to  edit  and  plot  contours  of 
various  components  of  stress  in  and  perpendicular  to  a given  plane. 


We  emphasize  that  the  key  to  the  development  of  such  a graphics 
package  is  in  fact  efficient  algorithms  for  the  numerical  inversion 
of  isoparametric  mappings. 
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2.  RESULTS  OF  INVESTIGATION 

A.  PRE-PROCESSOR:  VERIFICATION  OF  PROBLEM  GEOMETRY 

In  preparing  input  for  a multi-purpose  finite  element  program 
the  verification  that  the  three  dimensional  structure  has  been 
amenably  decomposed  into  "isoparametric''  curved  bricks  is  of  utmost 
importance.  As  a means  of  verifying  that  the  subdivision  contains 
no  anomalies,  it  is  desired  to  have  the  capability  of  passing  a 
given  plane  through  the  catenation  of  brick  elements  approximating 
the  structure  and  viewing  the  resulting  two-dimensional  intersection. 

Figure  2 illustrates  the  intersection  of  a three  dimensional 
isoparametric  brick  with  a plane.  This  figure  was  generated  using 
the  program  PLANIT,  C 11,123  developed  at  the  University  of  Pittsburgh. 
One  goal  of  this  work  program  was  the  development  of  computer  soft- 
ware to  plot  such  intersections  via  a CRT  and/or  an  incremental 
plotter  such  as  CALCOMP.  This  computer  program  determines  the 
intersection  of  a given  plane  with  a union  of  many  such  20-node 
isoparametric  elements.  Cf.  Figures  5 and  6. 

The  program  PLANIT  accepts  as  input  for  each  element 

(i)  element  number, 

(il)  coordinates  (xi,yi,z1)  of  the  nodes, 

(ill)  sufficient  information  to  determine  the  intersecting 
plane,  e.g.  point  and  normal. 

The  program  proceeds  to 

(!)  determine  those  elements  through  which  the  plane  passes; 

• » • • 

(ii)  define  the  curves  determined  by  the  intersection  of  the 
plane  with  each  element  from  (1); 


(ill)  produce  an  orthogonal  view  of  the  region  which  represents 
the  Intersection  of  the  plane  with  given  structure; 


FIGURE  2 


The  intersection  of  a plane  with  a single  20-node  brick 
Isometric  view,  (b)  View  along  normal  to  plane. 
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(iv)  store  the  information  necessary  to  retrieve  this  same 

plane  intersection  (if  desired)  to  be  used  in  the  post- 
processing mode  to  plot  various  stress  distributions. 


In  order  to  determine  the  Intersection  of  a plane 


(1)  a1x  + z2y  + a^z  + a^  « 0 

and  an  element  face  determined  by 


- “ 

- 

xu.n) 

* l (C,n) 

xi 

Y(5,n) 

yi 

ZU,n) 

± 

± 

_zi_ 

the  nonlinear  equation 

(3)  FU,n)  * a1X(5,n)  + a-YU,n)  + a^Z(C,n)  + * 0 

is  solved  numerically  for  n as  a function  of  5,  say  n * f(£)* 

Given  a specific  value  of  ( ■ (Q  then  F(5Q,n)  is  a quadratic  in  n 
which  is  easily  solved.  The  algorithm  tests  If  there  are  zero,  one, 
or  two  admissible  values  of  n.  Equation  (3)  is  obtained  by  sub- 
stituting (2)  into  (1).  The  locus  of  points  (X(5,f(5),  Y(£,f(C), 
Z(E»f(S))  determined  by  (3)  is  the  desired  curve  C in  3-space. 

If  n is  a unit  normal  to  the  given  plane  (Pand  P : (X  ,Y  ,Z  ) 

O 0 0 0 

is  a point  in  the  plane  designated  as  the  origin  then  a point 
is  also  in  the  plane  (P  if 

(P1-t50)  • n ■ 0. 

(As  input,  n and  a point  P1  in  the  plane  <P  could  be  given,  or, 
given  three  points  in  the  plane  n could  be'  computed.)  A coordinate 
system  is  established  by  the  orthogonal  unit  vectors 


el  ■ 


n x e. 


The  curve  C can  now  be  plotted  In  the  plane  2P  relative  to  the 


established  coordinate  system  by 


s - (P-PQ)  • e1 


t « (P-PQ)  • e2 


where  P traces  over  6. 


Figures  3 and  4 Illustrate  two  other  Intersections  of  single 
bricks  with  a plane.  As  indicated,  the  boundary  of  intersection 
may  be  disjoint  and  the  intersection  may  not  be  simply  connected, 
i.e.,  it  may  be  the  union  of  disjoint  lines  or  even  single  points. 
For  that  reason,  one  must  be  extremely  careful  in  connecting  the 
computed  points  of  the  boundary  of  the  intersection  by  a curve. 
Such  decisions  are  best  made  interactively  by  the  analyst  with 
the  aid  of  a graphics  terminal. 
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* The  pre-processor  strategy  Is  as  follows: 

• A geometry  file  for  the  proposed  finite  element  Idealization 
Is  constructed. 

• The  user  views  this  Idealization  and  chooses  a cross-sectional 
plane  of  the  structure  for  close  scrutiny. 

• The  user  specifies  all  or  some  of  those  brick  elements  he 
wishes  to  be  scanned. 

•The  Intersection  of  each  brick  with  the  specified  plane  Is 
determined  (if  such  exists)  and  displayed. 

•The  user  labels  each  portion  of  the  intersection  that  he  desires 
appended  to  a master  plot  file. 

•After  considering  the  totality  of  bricks  to  be  scrutinized, 
the  master  plot  file  is  displayed.  Voids  between,  or  overlapping 
of  elements  will  appear  as  unlabeled  or  multiple  labeled  subregions. 

•The  geometry  file  can  then  be  scanned  and  corrected  as  needed. 

Figure  5 illustrates  a portion  of  a reactor  vessel  near  the 
intersection  with  a 45°  lateral.  The  finite  element  idealization 
contains-  450  brick  elements.  This  figure  serves  to  demonstrate  the 
seemingly  Impossible  task  of  detecting  anomolles  In  the  geometric 
data  as  well  as  gauging  the  distribution  of  nodes.  The  above 
strategy  was  Implemented  for  this  structure  and  Figure  6 contains 
a plot  of  a plane  of  Intersection.  From  this  figure  It  Is  not 
difficult  to  see  that  an  error  exists  in  the  original  geometric 
data.  Specifically,  as  described  by  this  data,  elements  407  and 
408  are  not  contlnguous  to,  but  penetrate,  element  206. 


FIGURE  5.  Reactor  vessel  near  lateral'  nozzle 


FIGURE  6.  Plane  of  Intersection  and  a section 
of  that  Intersection  which  has  been  windowed 
for  greater  detail.  Note  the  overlap  between 
elements  206,  407  and  408.  There  are  two 
subregions  which  are  unlabeled  and  the  order  In 
which  their  boundaries  were  drawn  Indicates 
an  overlap. 


••• 
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B.  INVERSION  OF  TWO-DIMENSIONAL  ISOPARAMETRIC  MAPPINGS 

The  nature  of  two-dimensional  Isoparametric  mappings  ■* 

(where  J : (0,1)  * (0,1)  and  £ is  the  specified  element)  and  in 
particular  the  existence  of  inverses  for  such  mappings  has  been 
analyzed.  Various  results  are  reported  in  C83  and  include: 

(a)  The  straight  sided  quadralateral  mapping  T^  is  bijective 
(one-to-one)  if  and  only  if  £ is  convex.  (Cf.  Figure  7). 

(b)  The  8-node  quadratic  isoparametric  mapping  T2  is  bijec- 
tive if  its  Jacobian  is  nonzero  throughout  ,/  and  the  mapping  is 
one-to-one  on  the  boundary.  A class  of  element  shapes  is  specified 
for  which  the  Jacobian  is  nonzero.  (Cf.  Figure  8). 

(c)  For  the  8-node  quadratic  isoparametric  mapping  ?2» 
sufficient  conditions  are  given  to  guarantee  that  overspill  cannot 
occur.  Further,  under  rather  weak  assumptions  it  is  proven  that 
T2  is  bijective  if  and  only  if  overspill  does  not  occur. 

(d)  Perturbation  arguments  are  given  to  provide  another  set 
of  sufficient  conditions  for  the  bijectivity  of  Tj. 

Details  of  this  work  can  be  found  In  reference  [8], 

Mathematics  of  Computations.  32,  July  1978,  pp.  725-7^9. 


CONVEX  ELEMENT 
INVERTIBLE 


NCN-CONVEX  ELEMENT 
NON- INVERT IELE 


\ / 


FIGURE  7. 

Straight  aided  quadralateral 
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An  algorithm  for  numerically  Inverting  the  two-dimensional 
8-node  quadratic  isoparametric  mapping  has  been  developed.  The 
algorithm  is  based  on  bigradients  and  appears  to  be  relatively 
efficient.  (Cf.  Figure  9).  This  algorithm  will  determine  all 
pre-images  of  a point  in£,  including  those  which  do  not  lie  in 

A program  ( S0LV8 ) which  implements  this  algorithm  was  documented 
and  submitted  for  publication,  C103.  This  paper  was  revised  and 
resubmitted  in  June,  1979.  The  algorithm  now  includes  an  iterative 
improvement  strategy  based  on  Newton's  method. 


FIGURE  9 

The  pre-images  of  the  points  marked  "x"  were  found  to  be  those 
starked  "o". 
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C.  THE  CONSTANT  STRAIN  CONDITION 

The  constant  strain  condition  for  classes  of  second  and 
fourth  order  boundary  value  problems  was  investigated C Ph. D. 
Dissertation,  A.  E.  Prey,  19783.  Results  concerning  the  effect 
on  the  error  in  calculating  stiffness  matrices  when  using  elements 
with  curved  boundaries  include: 

(a)  Second  Order  Problems.  Exact  integration  preserves  an 
element’s  ability  to  satisfy  the  constant  strain  condition  even 
when  curved  boundaries  are  introduced.  For  any  reasonable  numeri- 
cal integration  formula,  it  is  shown  that  even  though  large  errors 
may  be  introduced  in  the  entries  of  the  stiffness  matrix  due  to 
rational  integrands  (caused  by  curved  element  boundaries)  the 
vectors  x and  j£  are  in  the  null  space  of  this  error  matrix. 

Hence  any  such  numerical  integration  formula  also  preserves  ar. 
element’s  ability  to.  reproduce  constant  strain.  This  work  is 
reported  in  C132,  which  has  been  accepted  for  publication. 

(b)  Fourth  Order  Problems.  We  were  unable  to  find  in  the 
literature  -a  ^-compatible  element  which  satisfies  the  constant 
strain  condition  for  elements  with  curved  boundaries.  Using  sub- 
parametric  element  mapping  concepts,  such  an  element  was  developed. 
(Paper  in  preparation.) 

D.  SUPPORT  PROGRAMS 

(1)  Computer  graphics  routines  BC8E2  and  B20E3  were  developed 
which  graphically  display  the  Images  of  two  and  three  dimensional 
isoparametric  mappings.  These  were  documented  in  C9J. 
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j*  (2)  An  algorithm  PLAINT,  [9]  was  developed  which  determines 

the  intersection  of  a single  20-node  isoparametric  (3-dimensional) 
brick  and  a given  plane.  The  necessary  software  for  graphic 
representation  of  this  intersection  was  also  developed.  (Cf. 

Figures  2,  3,  and  4.) 

(3)  A conversational  computer  graphics  program  PLANIT  was 
developed  using  PLAINT  to  plot  the  intersection  of  a given  plane 
with  a union  of  isoparametric  bricks.  Cf.  Figures  5 and  6. 

PLANIT  has  been  documented  in  [11],  and  was  reported  on  at  a 
symposium  [12]  and  in  Computers  and  Structures , 10,  1979,  pp.  149-154. 
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0 

0 

0 
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(4)  The  bigradient  algorithm  used  in  [10]  for  the  8-node 
element  was  considered  for  the  numerical  inversion  Of  the  three 
dimensional  20-node  isoparametric  element.  This  does  not  seem 
to  be  a fruitful  approach  due  to  the  complexity  of  the  elimination 
process.  An  algorithm  based  on  Newton's  method  was  implemented 
and  seems  to  perform  reasonably  well. 

E.  POST-PROCESSOR:  PLOTTING  OF  STRESS  DISTRIBUTIONS 

A finite  element  program  computes  various  components  of 
stress  in  terms  of  the  global  (x,y,z)  coordinate  system.  A 
program  STRSIT  [16]  was  developed  which  plots  in  a post-processing 
mode  contours  of  components  of  stress  in  and  perpendicular  to  a 
given  plane  which  intersects  the  structure. 

The  input  is  an  output  file  of  the  program  PLANIT.  This 
file  contains  the  intersection  of  the  plane  and  the  finite  element 


idealization.  In  addition,  the  element  geometry  file  is  available 


Aa  analyst  also  specifies  those  components  of  in-plane  and 


out-of-plane  stresses  which  he  desires  to  be  plotted. 

Given  the  mesh  point  displacements  (u. ,v  ,w. ) determined  by 


a finite  element  program,  the  displacement  components 


can  be  reconstructed  element-by-element.  The  displacement  vector 
T s (u,v,w)  is  independent  of  coordinate  system  and  can  be 
projected  onto  the  given  plane  and  a normal  to  the  plane  (cf . Figure 
10).  If  we  denote  the  new  components  of  6 relative  to  the  (r,s,t) 


FIGURE  10 


u,,v,,w',  then  the  strains  relative  to  the  (r,s,t)  directions  are 


These  computations  necessitate  inverting  various  isoparametric 
maps  and  their  associated  Jacobian  matrices.  Via  Hooke's  Law  the 


stresses  are  then  determined  as 
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where  D is  the  appropriate  6x6  elasticity  matrix. 

The  program  STRSIT  (element-by-element)  computes 
the  stress  state  tfrst  and  plots  those  components  desired  in  the 
plane  specified  in  the  input.  A sketch  of  such  a contour  plot  is 
given  in  Figure  11  where  the  plane  intersects  fourteen  different 
elements.  (The  program  HOLES4  [15 1 is  used  to  plot  the  contours). 


1 -0438787E-01 

2 0. 38034 1E+00 

3 0.804134E+00 

4 9.1227772*91 

9 0 163133E+01 

6 0.207499E-KJ1 

7 0 . 2496606+01 

8 0.2922212*91 


0 

0 

0 


FIGURE  11 

Contour  plot  of  a showing  also  the  Intersection  of  the 

s 

plane,  with  the  fourteen  bricks  which  contribute  to  the  stress 
distribution. 
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